function  [CS_i, CS_avg, E_i, E_avg, W_i, W_avg,  Error_i, Error_avg, ext_i, ext_avg, gvt_i, gvt_avg, SW_i, SW_avg, P_decision, P_decision_j] = Welfare_decomp_k_WTP_EEgap_parfor(param,param_FE,param_FEdemo,id_ij,J,p_elec,t_pigou,d_ext,elec,lifetime,r,MCPF,MT,price, price_denom, rebate_estar_denom, estar_denom, kwh, prod_denom,id_i, id_j, id_m, id_n, type_id2_3_demo1_denom, AV_demo1_denom,  ice_sc_demo1_denom, type_id2_3_demo2_denom,   AV_demo2_denom,  ice_sc_demo2_denom,type_id2_3_demo3_denom,   AV_demo3_denom,  ice_sc_demo3_denom,type_id2_3_demo4_denom,   AV_demo4_denom,  ice_sc_demo4_denom,type_id2_3_demo5_denom,  AV_demo5_denom,  ice_sc_demo5_denom, type_id2_3_demo6_denom,  AV_demo6_denom,  ice_sc_demo6_denom)
                                                                                                                                                

% Welfare formula:
% Delta CS = 1/eta *([log (sum_j U_j) +sum_j P_j(U^e-U)]_post_policy - [log (sum_j U_j) +sum_j P_j(U^e-U)]_pre_policy)

 %Parameters
  alpha=[0;param_FE(1:J-1,1)];
  %eta=param(1);
  eta=param(2);
  eta_experience=param(2);
  theta=param(3);
  tau=param(4);
  pi=param(5);
  beta=param_FEdemo;
  discount_fac=1/(1+r); 
  LLC=discount_fac*(1-discount_fac^lifetime)/(1-discount_fac); 
 
  FEprod_denom=alpha(prod_denom); 
 
 % kwh  
    kwh_denom=kwh(id_ij);
  
 % elec
    p_elec_ext=elec+p_elec+t_pigou;
    t_pigou_J=repmat(t_pigou,1,J);
  	elec_J=repmat(p_elec_ext,1,J);
    elec_J=elec_J';
    elec_denom=elec_J(id_ij).*kwh_denom;
    
       
 %Computations  
  				
 U_decision = FEprod_denom+eta*price_denom+eta*pi*rebate_estar_denom + eta*tau(1)*estar_denom...
                +eta*theta(1)*LLC*elec_denom...
  				+beta(1)*type_id2_3_demo1_denom +beta(2)*type_id2_3_demo2_denom + beta(3)*type_id2_3_demo3_denom...
  				+beta(4)*type_id2_3_demo4_denom +beta(5)*type_id2_3_demo5_denom + beta(6)*type_id2_3_demo6_denom...
				+beta(7)*AV_demo1_denom +beta(8)*AV_demo2_denom + beta(9)*AV_demo3_denom...
  				+beta(10)*AV_demo4_denom +beta(11)*AV_demo5_denom + beta(12)*AV_demo6_denom;				
 
 expU_decision=exp(U_decision);
  				
				
 U_experience= FEprod_denom+eta_experience*price_denom +eta_experience*LLC*elec_denom...
 				+eta_experience*tau*estar_denom + eta_experience*pi*rebate_estar_denom... 				
 				+beta(1)*type_id2_3_demo1_denom +beta(2)*type_id2_3_demo2_denom + beta(3)*type_id2_3_demo3_denom...
  				+beta(4)*type_id2_3_demo4_denom +beta(5)*type_id2_3_demo5_denom + beta(6)*type_id2_3_demo6_denom...
				+beta(7)*AV_demo1_denom +beta(8)*AV_demo2_denom + beta(9)*AV_demo3_denom...
  				+beta(10)*AV_demo4_denom +beta(11)*AV_demo5_denom + beta(12)*AV_demo6_denom;				
  					
  expU_experience=exp(U_experience);
    
  U_decision_ij = sparse(id_i,id_j,U_decision,id_m,id_n)';
  num_decision = sparse(id_i,id_j,expU_decision,id_m,id_n)';
  denom_decision=sum(num_decision,2);

  U_experience_ij = sparse(id_i,id_j,U_experience,id_m,id_n)';
  num_experience = sparse(id_i,id_j,expU_experience,id_m,id_n)';
  denom_experience=sum(num_experience,2);
  
  P_decision=num_decision./repmat(denom_decision,1,J);
  P_experience=num_experience./repmat(denom_experience,1,J);
  
  P_decision_j=sum(P_decision,1)/MT;
  P_experience_j=sum(P_experience,1)/MT;
  
  %Error_i=1000*abs(1/eta_experience)*sum(P_decision.*(U_experience_ij - U_decision_ij),2);
  %CS_i=1000*abs(1/eta_experience)*(log(denom_decision));
  Error_i=1000*abs(1/eta)*sum(P_decision.*(U_experience_ij - U_decision_ij),2);
  CS_i=1000*abs(1/eta)*(log(denom_decision));
  W_i=CS_i+Error_i;
  W_experience_i=1000*abs(1/eta_experience)*(log(denom_experience));
  
  Error_avg=sum(Error_i)/MT;
  CS_avg=sum(CS_i)/MT;
  W_avg=sum(W_i)/MT;
  W_experience_avg=sum(W_experience_i)/MT;
   
  E_ij=LLC*kwh*1000;
  E_i=sum(P_decision.*E_ij',2);
  E_avg=sum(E_i)/MT;
  E_experience_i=sum(P_experience.*E_ij',2);
  E_experience_avg=sum(E_experience_i)/MT;

  ext_ij=LLC*kwh.*d_ext*1000;
  ext_i=sum(P_decision.*ext_ij',2);
  ext_avg=sum(ext_i)/MT;
  ext_experience_i=sum(P_experience.*ext_ij',2);
  ext_experience_avg=sum(ext_experience_i)/MT;
  
  gvt_ij=LLC*kwh.*t_pigou_J'*1000;
  gvt_i=sum(P_decision.*gvt_ij',2);
  gvt_avg=sum(gvt_i)/MT;
  gvt_experience_i=sum(P_experience.*gvt_ij',2);
  gvt_experience_avg=sum(gvt_experience_i)/MT;
  
  SW_i=W_i-ext_i+MCPF*gvt_i;
  SW_avg=sum(SW_i)/MT;
  
  SW_experience_i=W_experience_i-ext_experience_i+MCPF*gvt_experience_i;
  SW_experience_avg=sum(SW_experience_i)/MT;
end 











